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Abstract 



O I A platelet collapse model of neutron starquakes is introduced. It 

-4— » ' displays self-organized criticality with a robust power-law behavior. 

C^ . The simulations indicate a near-constant exponent, whenever scaling 



is present. 
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Bak, Tang and Wiesenfeld [1] introduced the concept of self-organized 
criticality (SOC): large dissipative dynamical systems tend to drive them- 
selves to a critical state without the need to fine-tune system parameters. 
SOC seems to provide a natural explanation for the ubiquitous scale-free 
phenomena seen in nature. The critical state is a non-equilibrium state 
characterized by spatial correlation functions exhibiting power-law behavior. 
Until the finite size of the system comes into play, there is no intrinsic length 
scale in the dynamics, giving rise to the observed scaling. Thus the outstand- 
ing feature of SOC is the appearance of avalanches of dynamical activity of 
all sizes. Recently Bak [2] has given examples of SOC in astrophysics, includ- 
ing light emitted from quasars, emission of hard x-rays from solar flares and 
pulsar ghtches. In this paper we present a detailed model of pulsar glitches, 
which realizes SOC, characterized by partial or full collapse of the neutron 
star's surface induced by crustal strain. 

A pulsar glitch is an essentially instantaneous discontinuity of the rotating 
neutron star's period to shorter values. Morley and Garcia-Pelayo [3] have 
taken the world's pulsar glitch data and have shown empirically that the 
probability P{E) that a macro-glitch of energy E can occur in a glitching 
pulsar has the scaling law 

P{E) ~ E-^ (1) 

where the critical exponent, 7, due to the small number of recorded glitches, 
is not a well-determined number, lying in the range +1 < 7 < +2. 

There are many hundreds of known pulsars, of which only a minute num- 
ber actually glitch. Morley [4] has shown that a reasonable explanation of 
why pulsars glitch is the following: some fraction of pulsars are embedded 
in nebula or in regions of high interstellar mass density which allows them 
to accrete matter. All neutron stars and white dwarfs have an inherent in- 
stability induced by mass accretion, due to their inability to hydrostatically 
change their configuration. As one of these stars accretes matter and goes 
to a higher mass, its radius is also increasing which pushes the star to radial 
instabilityQ This mass accretion stress is added to the stress caused by the 
torque on the star, which causes the pulsar to be an ellipsode of revolution for 
the wrong angular speed. Together, these stresses induce tremendous strain 
in the neutron star's surface crust. When the threshold level is exceeded, the 



^Masses of neutron stars and white dwarfs in stable equilibrium are monotonic decreas- 
ing functions of their radius. 



crust relaxes, moving inward, resulting in a starquake. By conservation of 
angular momentum, this inward movement of surface crust causes the pulsar 
to immediately speed-up, the pulsar glitch. 

We now model this dynamical process. A pulsar has a misalignment 
of its spin axis S, with its magnetic dipole axis M (see fig. 1). In the M 
frame of reference, Ruderman [5] has argued that the surface of the star 
is divided into platelets which have varying surface sizes, constrained by 
the requirement that the total magnetic flux threaded through a platelet is 
approximately constant. We adopt the Olami et al. [6] cellular automaton 
model of earthquakes to this starquake model. Each platelet referenced (i,j) 
is connected to neighbors by mean shear moduli fis- The displacement of 
each platelet from its relaxed position is defined as dvij. The total radial 
stress on a given platelet is 

Sid = f^^s['2drij-dri^ij-dri+ij]+fi[''[2drij-drij^i-drij+i]+ficdrij . (2) 

In eq(2) /i'/, /i^" and ^c are respectively the longitudinal, latitudinal shear 
moduli and the compression modulus. The latter stress is associated with 
mass accretion, while the former are shear stresses. The total stress on a 
platelet increases uniformily (with a rate proportional to ^, where R is the 
star's radius) until one platelet reaches its threshold value and relaxes. This 
triggers the starquake. The redistribution of stress after a platelet collapse 
at position (i,j) is 

Si±l,j — ^ Si±l,j + ^Si±ij 

Si,j±i -^ Si,j±i + ^Sij±i (3) 



S,,, ^ 



where 



^^''^^' = 2/ii° + 2^/ + /./^'^- = ^"^^''^ ■ ^^^ 

For the purposes of this paper, we put /3'° = /?'" = p. Eqs(2)-(4) constitute 
the Olami et al. cellular automaton model [6]. We now add the important 
changes intrinsic to starquakes. We designate X{6i, (pj) as the surface platelet 



density function; that is, the platelet whose center is at position (i,j) has 
an area X{6i, (pj) in the M-system. Since the radial magnetic field is M/j = 
Mocos(^^) we parametrize \{9i,(f)j) as {9i, (pj are M system angles) 

The constant C is determined by the area normalization condition. We adopt 
a spherically shaped pulsar, and account for its eccentricity stresses directly 
as seen below. To handle the boundary topology of a sphere, we disconnect 
the pole and anti-pole points from their neighbors. By discretizing the 6j^ 
variable into M values, C satisfies CJ2ii2^ -^^i = 2(1 — Ax/ {a + 1)), where 

X = cos(6'jg), Axi = Ax = 2/M, Xi = ^^^r^- is discretized into A^ values 

by A0 = ^, so (pj = jA(f). There are thus A^ x (M - 2) platelets in the ill- 
system covering the star. In eq(5), a is the parameter which determines the 
ratio of the largest to smallest platelet sizes ~ ^r^- Since the eccentricities 
(e, eo) of pulsars are minute, we adopt a spherical configuration. To account 
for the centrifugal strain of having e deformation instead of eo (strain-free), we 
introduce a platelet threshold strain, Sth{ds)i which is a function of sin(^^), 
where dg is the polar angle in the S'-system: 

Sth{ds)i,j = Sril - a\ sm{9g)\ij) . (6) 

a in eq(6) is the parameter which accounts for the reduced threshold for 
platelet collapse induced by the centrifugal strain. All stresses are in units 
of St, which can be put equal to 1. By geometry | sin(6'^)|jj = (1 — 



(sin(-?/') cos(0j) Y 1 — xf + Xi cos(-?/'))^)^/^, where tp is the angle between S and 

M. The last major change is the new boundary conditions. All platelets 
except those ringing the pole and anti-pole have continuous boundary con- 
ditions of four neighbors. Those ringing the poles have only three. The 
parameters of the model are ip,a,a,j3 and the number A^ x (M — 2). The 
simulation proceeds as follows: 



1) Initalize all platelet stresses to a random value between and S*, 



thyvgjij. 



2) If any Sij > Sth{Og)ij, the platelet relaxes. Redistribute the stress on 
Sij according to eq(3) and add the areas of the collapsed platelets. 



3) Repeat step 2) until the starquake is fully evolved. 

4) Locate the platelet with the largest strain such that AS = Sth{Og)ij — 
Si J is smallest. 

5) Add AS* to all platelets (uniform mass accretion) and return to step 
2). 

The energy of a starquake comes from gravity. We assume the platelets fall 
a fixed relative distance -^, so the energy is just proportional to the area 
of the collapsed platelets: E ~ J2coiiapsedi,j ^{^iy(pj)- Interestingly enough, 
there are rare starquakes where more than the total surface area collapses 
due to the near-threshold level of all platelets and the continuous boundary 
conditions. 

In fig. 2, we present log ^ versus logE, with the listed parameters^. 
No SOC is present. This can be understood rather simply: the value of a 
means there is a disproportionate area size distribution among platelets which 
means that the starquakes are dominated by the large platelet collapses. By 
changing a to a value where the area density is not greatly varying, fig. 3, 
SOC is realized. Simulations indicate that SOC is present for all small a 
values and the exponent is almost independent of the values of ip and a. 
Only variation in (3 produces a weak modification of the scaling exponent, 
fig. 4. 

Our computer results strongly support SOC in the phenomena of pulsar 
glitches. They indicate that mass accretion is the critical element driving the 
neutron star instability. The critical exponent should lie between 1.5 < 7 < 
2, depending only on the ratios of shear and compression moduli. 



^AU runs had 999999 starquakes. 
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Figure Captions 

1. Platelets on the rotating neutron star with misaligned spin S and mag- 
netic dipole M axis. 

2. log^ versus logE, for ■^=0.5, a=0.01, a=.l, N = M = 50 and /5=0.2. 

3. Same as fig. 2, except a=2.0. 

4. Critical exponent 7 versus stress parameter (3. 
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